home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / eigen_ii.pro < prev    next >
Text File  |  1997-07-08  |  5KB  |  151 lines
  1. ; $Id: eigen_ii.pro,v 1.5 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ; Copyright (c) 1993-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5.  
  6. ;+
  7. ; NAME:
  8. ;       EIGEN_II
  9. ;
  10. ; PURPOSE:
  11. ;       This function computes the eigenvectors of an N by N 
  12. ;       real, nonsymmetric matrix. The result is an array of 
  13. ;       complex type with a column dimension equal to the 
  14. ;       number of eigenvalues and a row dimension equal to N.
  15. ;
  16. ; CATEGORY:
  17. ;       Numerical Linear Algebra.
  18. ;
  19. ; CALLING SEQUENCE:
  20. ;       Result = EIGEN_II(A, Eval)
  21. ;
  22. ; INPUTS:
  23. ;           A:  An N by N matrix real, nonsymmetric matrix.
  24. ;        Eval:  An n-element complex vector of eigenvalues.    
  25. ;
  26. ; KEYWORD PARAMETERS:
  27. ;      Double:  If set to a non-zero value, computations are done in
  28. ;               double precision arithmetic.
  29. ;               NOTE: Since IDL lacks a double-precision complex 
  30. ;                     data type, computations are done internally
  31. ;                     in double-precision and the result is truncated
  32. ;                     to single-precision complex.
  33. ;       Itmax:  The number of iterations performed in the computation
  34. ;               of each eigenvector. The default value is 4.
  35. ;
  36. ; EXAMPLE:
  37. ; 1) A real, nonsymmetric matrix with real eigenvalues/eigenvectors.
  38. ;       Define an N by N real, nonsymmetric array.
  39. ;         a = [[  7.3, 0.2, -3.7], $
  40. ;              [-11.5, 1.0,  5.5], $
  41. ;              [ 17.7, 1.8, -9.3]]
  42. ;       Transpose the array into IDL column (matrix) format.
  43. ;         at = transpose(a)
  44. ;       Compute the eigenvalues of a.
  45. ;         eval = NR_HQR(NR_ELMHES(at))
  46. ;       Print the eigenvalues.
  47. ;         print, eval
  48. ;       Compute the eigenvectors of a.
  49. ;         evec = EIGEN_II(at, eval)
  50. ;       Print the eigenvectors.
  51. ;         print, evec(0,*), evec(1,*), evec(2,*)
  52. ;       Check the accuracy of each eigenvalue/eigenvector (lamda/x) 
  53. ;       pair using the mathematical definition of an eigenvector;
  54. ;       Ax - (lambda)x = 0
  55. ;         print, (evec(0,*) # a) - (eval(0) * evec(0,*))
  56. ;         print, (evec(1,*) # a) - (eval(1) * evec(1,*))
  57. ;         print, (evec(2,*) # a) - (eval(2) * evec(2,*)) 
  58. ; 2) A real, nonsymmetric matrix with complex eigenvalues/eigenvectors.
  59. ;       Define an N by N real, nonsymmetric array.
  60. ;       a = [[ 0.0, 1.0], $
  61. ;            [-1.0, 0.0]]
  62. ;       Transpose the array into IDL column (matrix) format.
  63. ;         at = transpose(a)
  64. ;       Compute the eigenvalues of a.
  65. ;         eval = NR_HQR(NR_ELMHES(at))
  66. ;       Print the eigenvalues.
  67. ;         print, eval
  68. ;       Compute the eigenvectors of a.
  69. ;         evec = EIGEN_II(at, eval, /double)
  70. ;       Print the eigenvectors.
  71. ;         print, evec(0,*), evec(1,*)
  72. ;       Check the accuracy of each eigenvalue/eigenvector (lamda/x)
  73. ;       pair using the mathematical definition of an eigenvector;
  74. ;       Ax - (lambda)x = 0
  75. ;         print, (evec(0,*) # a) - (eval(0) * evec(0,*))
  76. ;         print, (evec(1,*) # a) - (eval(1) * evec(1,*))
  77. ;
  78. ; PROCEDURE:
  79. ;       EIGEN_II.PRO computes the set of eigenvectors that correspond
  80. ;       to a given set of eigenvalues using Inverse Subspace Iteration.
  81. ;       The eigenvectors are computed up to a scale factor and are of
  82. ;       Euclidean length.
  83. ;       The functions NR_ELMHES and NR_HQR may be used to find the  
  84. ;       eigenvalues of an N by N matrix real, nonsymmetric matrix.
  85. ;
  86. ; REFERENCE:
  87. ;       Numerical Recipes, The Art of Scientific Computing (Second Edition)
  88. ;       Cambridge University Press
  89. ;       ISBN 0-521-43108-5
  90. ;
  91. ; MODIFICATION HISTORY:
  92. ;           Written by:  GGS, RSI, December 1993
  93. ;-
  94.  
  95. function eigen_ii, a, eval, double = double, itmax = itmax
  96.  
  97.   on_error, 2  ;return to caller if error occurs.
  98.  
  99.   dim = size(a)
  100.   if dim(1) ne dim(2) then stop, $
  101.     ' EIGEN_II: A is not square.'
  102.  
  103. ; Set default values for keyword parameters.
  104.   if keyword_set(double) eq 0 then double =  0
  105.   if keyword_set(itmax)  eq 0 then itmax  =  4
  106.  
  107.   enum = n_elements(eval)            ;Number of eigenvalues.
  108.   evec = complexarr(enum, dim(1))    ;Eigenvector storage array.
  109.   diag = indgen(dim(1)) * (dim(1)+1) ;Diagonal indices of a. 
  110.  
  111.   for k = 0, enum - 1 do begin
  112.     alud = a  ;Create a new copy of alud for next eigenvalue computation.
  113.  
  114. ;   Complex eigenvalue ?
  115.     if imaginary(eval(k)) ne 0 then begin
  116.       alud = complex(alud)
  117.       alud(diag) = alud(diag) - eval(k)
  118.       re = float(alud)
  119.       im = imaginary(alud)
  120.       comp = [[ re, im], $
  121.               [-im, re]]
  122.       b = replicate(1.0, 2*dim(1))
  123.       b = b / sqrt(total(b^2, 1))
  124.       ludc, comp, index, double = double, /column
  125.       it = 0
  126.       while it lt itmax do begin
  127.         x = lusol(comp, index, b, double = double, /column)
  128.         b = x / sqrt(total(x^2, 1))  ;Normalize eigenvector.
  129.         it = it + 1
  130.       endwhile
  131.       evec(k, *) = complex(b(0:dim(1)-1), b(dim(1):*))
  132.  
  133.     endif else begin
  134. ;   Real eigenvalue !
  135.       alud(diag) = alud(diag) - float(eval(k))
  136.       b = replicate(1.0, dim(1))
  137.       b = b / sqrt(total(b^2, 1))
  138.       ludc, alud, index, double = double, /column
  139.       it = 0
  140.       while it lt itmax do begin
  141.         x = lusol(alud, index, b, double = double, /column)
  142.         b = x / sqrt(total(x^2, 1))  ;Normalize eigenvector.
  143.         it = it + 1
  144.       endwhile
  145.       evec(k, *) = complex(b, 0.0)
  146.     endelse
  147.   endfor
  148.   return, evec
  149. end
  150.  
  151.